# Author: Mark Richardson
# Purpose: Execute analysis in the appendix

# Load packages
library(dplyr)
library(kableExtra)
library(fixest)
library(ggplot2)
library(margins)
library(patchwork)
library(stringr)

#### Table A.1 ####

joint <- read.table(file = "data/table_a_1.csv",
                    header = TRUE)

kbl(joint %>% rename(`Rarely-Never` = Rarely.Never), booktabs = TRUE, linesep = "",
    caption = "Joint Distribution of Strength and Frequency", escape = FALSE,
    digits = 2) %>%
  kable_styling() %>%
  footnote(general = "Note: N = 1,368; Cell entries are proportions.", threeparttable = TRUE, footnote_as_chunk = TRUE, general_title = "")



#### Figure A.1 ####

dis_means <- readr::read_csv(file = "data/fig_a_1_dis_means.csv")

ggplot(dis_means, aes(pty_dis_freq_jitter, pty_dis_strg_jitter)) +
  geom_point(alpha = 1/2) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_x_continuous(breaks = 1:4, limits = c(0.9, 4.1),
                     labels = c("Never/Rarely", "Sometimes", "Often", "Always")) +
  scale_y_continuous(breaks = 1:4, limits = c(0.9, 4.1),
                     labels = c("No disagreement", "Low", "Moderate", "High")) +
  labs(x = "Frequency of Disagreement",
       y = "Intensity of Disagreement",
       caption = "Plot inlcudes a fitted line from a bivariate OLS model.\nPoints are jittered to improve readiblity and prevent inference of individual-level responses.") +
  theme_bw()


# clean up
rm(dis_means, joint)

#### Figures A.2, A.3 and A.4 ####

# Load data
load("data/02_mis_data.RData")
load("data/05_dis_data_w_traits.RData")

# Preliminary wrangling
mis_data <- mis_data %>%
  arrange(pty_dis) %>%
  mutate(plot_index = 1:n())

dis_data_wt <- dis_data_wt %>%
  mutate(ideo_cat = case_when(ideo_rating < 0 & ideo_ub < 0 ~ "Liberal",
                              ideo_rating > 0 & ideo_lb > 0 ~ "Conservative",
                              TRUE ~ "Moderate")) %>%
  mutate(ideo_cat = if_else(is.na(ideo_rating), "No measure", ideo_cat),
         ideo_cat = factor(ideo_cat, levels = c("Conservative", "Moderate", "Liberal", "No measure"))) %>%
  arrange(pty_dis) %>%
  mutate(plot_index = 1:n())


dis_data_wt <- dis_data_wt %>%
  mutate(agency_plot = str_replace_all(agency, "United States", "U.S."),
         agency_plot = if_else(!is.na(dept_acr), paste0(agency_plot, " (", dept_acr, ")"), agency_plot),
         agency_plot = str_replace_all(agency_plot, "Office of the Secretary", "Office of the Sec."))




# Figure A.2

ggplot(dis_data_wt) +
  geom_point(aes(pty_dis, plot_index, shape = ideo_cat, color = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -5.5, y = plot_index,
                label = agency_plot,
                hjust = 1)) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12)) +
  labs(x = "Latent Partisan Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:4, limits = c(-17, 6), sec.axis = dup_axis()) +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  geom_vline(xintercept = 0, linetype = "dashed")

# Figure A.3

mis_data <- mis_data %>%
  arrange(pty_dis) %>%
  mutate(plot_index = 1:n())

ggplot(mis_data) +
  geom_point(aes(pty_dis, plot_index)) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index)) +
  geom_text(aes(x = -1, y = plot_index,
                label = str_replace(str_replace(mission, "Services", "Svcs."), "Education", "Educ."),
                hjust = 1)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "Latent Disagreement",
       y = "",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -1:1, limits = c(-3.5, 1), sec.axis = dup_axis()) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12))


# Differences in mission-level latent disagreement (pg. 12 in the appendix)

load("data/04_mis_draws.RData")

health_v_transpo <- (draws_mu %>% pull(`mu[9]`)) - draws_mu %>% pull(`mu[14]`)

hvt_bounds <- quantile(health_v_transpo, probs = c(0.025, 0.50, 0.975))

env_v_def<- (draws_mu %>% pull(`mu[13]`)) - draws_mu %>% pull(`mu[12]`)

evd_bounds <- quantile(env_v_def, probs = c(0.025, 0.50, 0.975))

# Health v. Transportatoin
mean(health_v_transpo) %>% format(digits = 2, nsmall = 2)

hvt_bounds[1]  %>% format(digits = 2, nsmall = 2)

hvt_bounds[3]  %>% format(digits = 2, nsmall = 2)

# Natural Resources v. Environment and National Defense
mean(env_v_def) %>% format(digits = 2, nsmall = 2)

evd_bounds[1]  %>% format(digits = 2, nsmall = 2)

evd_bounds[3]  %>% format(digits = 2, nsmall = 2)

# Figure A.4 preliminaries

msns <- list()

for (i in 1:nrow(mis_data)) {
  
  msns[[i]] <- bind_rows(mis_data %>%
                           slice(i) %>%
                           mutate(ideo_cat = "No measure") %>%
                           rename(agency_plot = mission),
                         dis_data_wt %>% filter(mission == mis_data$mission[i]) %>%
  mutate(agency_plot = str_replace_all(agency, "United States", "U.S."),
         agency_plot = if_else(!is.na(dept_acr), paste0(agency_plot, " (", dept_acr, ")"), agency_plot),
         agency_plot = str_replace_all(agency_plot, "Office of the Secretary", "Office of the Sec.")))
  
    msns[[i]] <- msns[[i]] %>%
    arrange(pty_dis) %>%
    mutate(plot_index = 1:n())
  
}

x_axis_lb <- -9


# Figure A.4 plots

i <- 1

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 2

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 3

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 4

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 5

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 6

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 7

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 8

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "bottom") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 9

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 10

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 11

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 12

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 13

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 14

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "none") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


i <- 15

ggplot(msns[[i]]) +
  geom_point(aes(pty_dis, plot_index, color = ideo_cat, shape = ideo_cat), size = 2, alpha = 1) +
  geom_segment(aes(x = pty_dis_lb_95, xend = pty_dis_ub_95,
                   y = plot_index, yend = plot_index,
                   color = ideo_cat), alpha = 0.5) +
  geom_text(aes(x = -4, y = plot_index,
                label = agency_plot),
                hjust = 1, size = 4) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(size = 12),
        legend.position = "bottom") +
  labs(title = paste0("Mission Area: ", na.omit(msns[[i]]$mission)),
       x = "Latent Disagreement",
       y = "",
       color = "Agency Ideology",
       shape = "Agency Ideology",
       caption = "Horizontal lines are the 95% Bayesian credible interval.") +
  scale_x_continuous(breaks = -4:3, limits = c(x_axis_lb , 3), sec.axis = dup_axis()) +
  scale_color_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_shape_manual(limits = c("Conservative", "Moderate", "Liberal", "No measure"),
                     values = 15:18) +
  geom_vline(xintercept = 0, linetype = "dashed")


# clean up
rm(draws_mu, mis_data, msns, env_v_def, evd_bounds, health_v_transpo, hvt_bounds,
   i, x_axis_lb)

#### Section A.6.1 and Figure A.5 ####

# Load agenda_ch_means
load("data/07_agenda_ch_means.RData")

dis_ch <- full_join(dis_data_wt, agenda_ch_out,
                    by = c("agency", "agency_index"))

ggplot(dis_ch) +
  geom_smooth(aes(pty_dis, ch_agenda_jitter), method = "lm") +
  geom_point(aes(pty_dis, ch_agenda_jitter, shape = ideo_cat, color = ideo_cat), size = 2, alpha = 1) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "EPA"), data = dis_ch %>% filter(agency == "Environmental Protection Agency"),
            nudge_x = 0.1, nudge_y = 0.15) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "BFS"), data = dis_ch %>% filter(agency == "Bureau of the Fiscal Service"),
            nudge_x = 0.25) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "Army &\nNavy"), data = dis_ch %>% filter(agency == "Army"),
            nudge_x = 1, nudge_y = -0.9) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "USAF"), data = dis_ch %>% filter(agency == "Air Force"),
            nudge_x = -0.3, nudge_y = 0) +
  geom_segment(aes(x = 1.1, y = 2.2, xend = 0.3, yend = 2.75), color = "gray", lineend = "round", arrow = arrow(length = unit(0.1, "inches"))) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "BLS"), data = dis_ch %>% filter(agency == "Bureau of Labor Statistics"),
            nudge_x = 0, nudge_y = -0.1) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "NTSB"), data = dis_ch %>% filter(agency == "National Transportation Safety Board"),
            nudge_x = 0, nudge_y = 0.1) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "USCIS"), data = dis_ch %>% filter(agency == "United States Citizenship and Immigration Services"),
            nudge_x = -0.2, nudge_y = 0.15) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "USCBP"), data = dis_ch %>% filter(agency == "United States Customs and Border Protection"),
            nudge_x = 0, nudge_y = -0.1) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "ATF"), data = dis_ch %>% filter(agency == "Bureau of Alcohol, Tobacco, Firearms and Explosives"),
            nudge_x = 0, nudge_y = 0.1) +
  geom_text(aes(pty_dis, ch_agenda_jitter, label = "NIST"), data = dis_ch %>% filter(agency == "National Institute of Standards and Technology"),
            nudge_x = 0, nudge_y = -0.1) +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_y_continuous(breaks = 1:4, labels = c("None", "Minimal", "Moderate", "Significant"), limits = c(1, 4.25)) +
  labs(x = "Latent Partisan Disagreement", y = "Mean Magnitude of Agenda Change",
       shape = "Agency Ideology", color = "Agency Ideology",
       caption = "Includes a fitted line from a bivariate OLS regression.\nShaded area gives the 95% confidence interval.\nPoints are jittered to improve readability and prevent inference of\nindividual-level responses.") +
  theme_bw()

# R^2
summary(lm(ch_agenda_jitter ~ pty_dis, data = dis_ch))$r.squared

# clean up
rm(dis_ch, agenda_ch_out)

#### Section A.6.2 and Figure A.6 ####

ggplot(dis_data_wt) +
  geom_smooth(aes(pty_dis, skills_rating), method = "lm") +
  geom_point(aes(pty_dis, skills_rating, shape = ideo_cat, color = ideo_cat), size = 2, alpha = 1) +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  labs(x = "Latent Partisan Disagreement", y = "Workforce Skill",
       shape = "Agency Ideology", color = "Agency Ideology",
       caption = "Includes a fitted line from a bivariate OLS regression.\nShaded area gives the 95% confidence interval.") +
  theme_bw()

# Model discussed in Section A.6.2
sk <- feols(skills_rating ~ pty_dis,
            data = dis_data_wt,
            vcov = "hetero")
 
etable(sk)
nobs(sk)
 
# Clean up
rm(sk)

#### Section A.6.3 ####

load("data/06_confirmation_votes.RData")

noms <- noms %>%
  mutate(abs_ideo = abs(ideo_rating))

mn <- noms %>%
  group_by(agency) %>%
  summarize(n_votes = sum(!is.na(pty_dif_vv)),
            pty_dis = mean(pty_dis),
            ideo_rating = mean(ideo_rating),
            pty_dif_vv = mean(pty_dif_vv, na.rm = TRUE),
            abs_ideo = mean(abs_ideo),
            ideo_cat = first(ideo_cat))



# Discussion of data on page 22

# Prop dems
table(noms$yea_dem/(noms$yea_dem + noms$nay_dem))

# Dems voting nay
table(noms$nay_dem)

# Prop reps
table(noms$yea_rep/(noms$yea_rep + noms$nay_rep))

# Number of votes (including voice votes)
table(is.na(noms$pty_dif_vv))

# Number of agencies with a vote
noms %>% filter(!is.na(pty_dif_vv)) %>% distinct(agency) %>% nrow()

# Votes per agency
summary(mn$n_votes)

# Votes per agency
mn %>% filter(agency == "Army") %>% pull(n_votes)
mn %>% filter(agency == "Air Force") %>% pull(n_votes)
mn %>% filter(agency == "Navy") %>% pull(n_votes)
mn %>% filter(agency == "United States Customs and Border Protection") %>% pull(n_votes)
mn %>% filter(agency == "Bureau of Alcohol, Tobacco, Firearms and Explosives") %>% pull(n_votes)

#### Confirmation models ####

# Create variable names for tables
setFixest_dict(pty_dis = "Partisan Dis.",
               ideo_rating = "Agency Ideology",
               agency = "Agency",
               pty_dif_vv = "Difference in Party Support",
               voice_vote = "Voice Vote",
               pty_unity_vv = "Party Unity Vote",
               abs_ideo = "|Agency Ideology|")


# OLS models of partisan difference in support (vote-level) (Table A.2)
lm_vv <- feols(pty_dif_vv ~ pty_dis,
               data = noms,
               cluster = ~agency)

lm_id <- feols(pty_dif_vv ~ ideo_rating,
               data = noms,
               cluster = ~agency)

lm_id_abs <- feols(pty_dif_vv ~ abs_ideo,
                   data = noms,
                   cluster = ~agency)

lm_vv_id <- feols(pty_dif_vv ~ pty_dis + ideo_rating,
                  data = noms,
                  cluster = ~agency)

lm_vv_id_int <- feols(pty_dif_vv ~ pty_dis * ideo_rating,
                      data = noms,
                      cluster = ~agency)

lm_vv_id_int_abs <- feols(pty_dif_vv ~ pty_dis * abs_ideo,
                          data = noms,
                          cluster = ~agency)

# OLS models of partisan difference in support (agency-level) (Table A.3)
lm_vv_mn <- feols(pty_dif_vv ~ pty_dis,
                  data = mn,
                  vcov = "hetero")

lm_id_mn <- feols(pty_dif_vv ~ ideo_rating,
                  data = mn,
                  vcov = "hetero")

lm_id_abs_mn <- feols(pty_dif_vv ~ abs_ideo,
                      data = mn,
                      vcov = "hetero")

lm_vv_id_mn <- feols(pty_dif_vv ~ pty_dis + ideo_rating,
                     data = mn,
                     vcov = "hetero")

lm_vv_id_int_mn <- feols(pty_dif_vv ~ pty_dis * ideo_rating,
                         data = mn,
                         vcov = "hetero")

lm_vv_id_int_abs_mn <- feols(pty_dif_vv ~ pty_dis * abs_ideo,
                             data = mn,
                             vcov = "hetero")


# Prob. of party unity vote (Table A.4)

# Models for etable

pr_uty_dis <- feglm(pty_unity_vv ~ pty_dis,
                data = noms,
                family = binomial(link = "probit"),
                cluster = ~agency)

pr_uty_id <- feglm(pty_unity_vv ~ ideo_rating,
                   data = noms,
                   family = binomial(link = "probit"),
                   cluster = ~agency)

pr_uty_id_abs <- feglm(pty_unity_vv ~ abs_ideo,
                   data = noms,
                   family = binomial(link = "probit"),
                   cluster = ~agency)

pr_uty_dis_id <- feglm(pty_unity_vv ~ pty_dis + ideo_rating,
                       data = noms,
                       family = binomial(link = "probit"),
                       cluster = ~agency)

pr_uty_dis_id_int <- feglm(pty_unity_vv ~ pty_dis * ideo_rating,
                           data = noms,
                           family = binomial(link = "probit"),
                           cluster = ~agency)

pr_uty_dis_id_int_abs <- feglm(pty_unity_vv ~ pty_dis * abs_ideo,
                           data = noms,
                           family = binomial(link = "probit"),
                           cluster = ~agency)


# Models for margins (margins package does not accept fixest objects)

pr_uty_dis_mg <- glm(pty_unity_vv ~ pty_dis,
                data = noms,
                family = binomial(link = "probit"))

pr_uty_id_mg <- glm(pty_unity_vv ~ ideo_rating,
                   data = noms,
                   family = binomial(link = "probit"))

pr_uty_id_abs_mg <- glm(pty_unity_vv ~ abs_ideo,
                   data = noms,
                   family = binomial(link = "probit"))

pr_uty_dis_id_mg <- glm(pty_unity_vv ~ pty_dis + ideo_rating,
                       data = noms,
                       family = binomial(link = "probit"))

pr_uty_dis_id_int_mg <- glm(pty_unity_vv ~ pty_dis * ideo_rating,
                           data = noms,
                           family = binomial(link = "probit"))

pr_uty_dis_id_int_abs_mg <- glm(pty_unity_vv ~ pty_dis * abs_ideo,
                           data = noms,
                           family = binomial(link = "probit"))




# Prob. of voice vote (Table A.5)

# Models for etable

pr_vv_dis <- feglm(voice_vote ~ pty_dis,
                data = noms,
                family = binomial(link = "probit"),
                cluster = ~agency)

pr_vv_id <- feglm(voice_vote ~ ideo_rating,
                   data = noms,
                   family = binomial(link = "probit"),
                   cluster = ~agency)

pr_vv_id_abs <- feglm(voice_vote ~ abs_ideo,
                   data = noms,
                   family = binomial(link = "probit"),
                   cluster = ~agency)

pr_vv_dis_id <- feglm(voice_vote ~ pty_dis + ideo_rating,
                       data = noms,
                       family = binomial(link = "probit"),
                       cluster = ~agency)

pr_vv_dis_id_int <- feglm(voice_vote ~ pty_dis * ideo_rating,
                           data = noms,
                           family = binomial(link = "probit"),
                           cluster = ~agency)

pr_vv_dis_id_int_abs <- feglm(voice_vote ~ pty_dis * abs_ideo,
                           data = noms,
                           family = binomial(link = "probit"),
                           cluster = ~agency)


# Models for margins

pr_vv_dis_mg <- glm(voice_vote ~ pty_dis,
                data = noms,
                family = binomial(link = "probit"))

pr_vv_id_mg <- glm(voice_vote ~ ideo_rating,
                   data = noms,
                   family = binomial(link = "probit"))

pr_vv_id_abs_mg <- glm(voice_vote ~ abs_ideo,
                   data = noms,
                   family = binomial(link = "probit"))

pr_vv_dis_id_mg <- glm(voice_vote ~ pty_dis + ideo_rating,
                       data = noms,
                       family = binomial(link = "probit"))

pr_vv_dis_id_int_mg <- glm(voice_vote ~ pty_dis * ideo_rating,
                           data = noms,
                           family = binomial(link = "probit"))

pr_vv_dis_id_int_abs_mg <- glm(voice_vote ~ pty_dis * abs_ideo,
                           data = noms,
                           family = binomial(link = "probit"))




#### Figure A.7 ####

id <- ggplot(noms, aes(ideo_rating, pty_dif_vv)) +
  geom_point(aes(color = ideo_cat, shape = ideo_cat)) +
  geom_smooth(method = "lm") +
  labs(x = "Agency Ideology",
       y = "Difference in Party Support",
       color = "Agency Ideology",
       shape = "Agency Ideology") +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits =  c(-0.3, 1)) +
  theme_bw() +
  theme(legend.position = "none")

dis <- ggplot(noms, aes(pty_dis, pty_dif_vv)) +
  geom_point(aes(color = ideo_cat, shape = ideo_cat)) +
  geom_smooth(method = "lm") +
  labs(title = "Vote-level",
       x = "Latent Partisan Disagreement",
       y = "Difference in Party Support",
       color = "Agency Ideology",
       shape = "Agency Ideology") +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits =  c(-0.3, 1)) +
  theme_bw() +
  theme(legend.position = "none")

id_mn <- ggplot(mn, aes(ideo_rating, pty_dif_vv)) +
  geom_point(aes(color = ideo_cat, shape = ideo_cat)) +
  geom_smooth(method = "lm") +
  labs(x = "Agency Ideology",
       y = "Difference in Party Support",
       color = "Agency Ideology",
       shape = "Agency Ideology") +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits =  c(-0.3, 1)) +
  theme_bw() +
  theme(legend.position = "none")

dis_mn <- ggplot(mn, aes(pty_dis, pty_dif_vv)) +
  geom_point(aes(color = ideo_cat, shape = ideo_cat)) +
  geom_smooth(method = "lm") +
  labs(title = "Agency-level",
       x = "Latent Partisan Disagreement",
       y = "Difference in Party Support",
       color = "Agency Ideology",
       shape = "Agency Ideology") +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), limits =  c(-0.3, 1)) +
  theme_bw()

(dis + dis_mn) / (id + id_mn) +
  plot_annotation(caption = "Includes a fitted line from a bivariate OLS regression. Shaded area gives the 95% confidence interval.") +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")
  



#### Table A.2 ####

etable(lm_vv, lm_id, lm_id_abs, lm_vv_id, lm_vv_id_int, lm_vv_id_int_abs,
       tex = TRUE, digits = "r2", digits.stats = "r2",
       order = "!Constant", fitstat = c("n", "r2", "ar2"),
       style.tex = style.tex(stats.title = "\\midrule"),
       title = "OLS Models of Difference in Party Support (Vote-level)",
       label = "pty_dif_table_vv")


#### Table A.3 ####

etable(lm_vv_mn, lm_id_mn, lm_id_abs_mn, lm_vv_id_mn, lm_vv_id_int_mn,
       lm_vv_id_int_abs_mn,
       tex = TRUE, digits = "r2", digits.stats = "r2",
       order = "!Constant", fitstat = c("n", "r2", "ar2"),
       style.tex = style.tex(stats.title = "\\midrule"),
       title = "OLS Models of Difference in Party Support (Agency-level)",
       label = "pty_dif_table_mn")



# Confidence intervals models 1 and 2 in Table A.2

confint(lm_vv)

confint(lm_id)

# Army v. EPA difference in partisan disagreement
round(mn %>% filter(agency == "Environmental Protection Agency") %>% pull(pty_dis) -
        mn %>% filter(agency == "Army") %>% pull(pty_dis), digits = 2)

# Army v. EPA difference in agency ideology
round(mn %>% filter(agency == "Army") %>% pull(ideo_rating) -
        mn %>% filter(agency == "Environmental Protection Agency") %>% pull(ideo_rating), digits = 2)

# partisan disagreement calc
1.75 * 0.24

# ideology calc
3.45 * -0.08
 


#### Table A.4 ####

etable(pr_uty_dis, pr_uty_id, pr_uty_id_abs,
       pr_uty_dis_id, pr_uty_dis_id_int, pr_uty_dis_id_int_abs,
       tex = TRUE, digits = "r2", digits.stats = "r2",
       order = "!Constant", fitstat = c("n", "pr2", "aic"),
       style.tex = style.tex(stats.title = "\\midrule"),
       title = "Probit Models of Party Unity Votes")



#### Table A.5 ####

etable(pr_vv_dis, pr_vv_id, pr_vv_id_abs,
       pr_vv_dis_id, pr_vv_dis_id_int, pr_vv_dis_id_int_abs,
       tex = TRUE, digits = "r2", digits.stats = "r2",
       order = "!Constant", fitstat = c("n", "pr2", "aic"),
       style.tex = style.tex(stats.title = "\\midrule"),
       title = "Probit Models of Voice Votes")



#### Predicted probabilities for Figure A.8 ####

pr_unity <- summary(prediction(pr_uty_dis_mg,
                    vcov = vcov_cluster(pr_uty_dis, cluster = ~agency),
                    at = list(pty_dis = unique(pr_uty_dis_mg$model$pty_dis))))

pr_vv <- summary(prediction(pr_vv_dis_mg,
                 vcov = vcov_cluster(pr_vv_dis, cluster = ~agency),
                 at = list(pty_dis = unique(pr_vv_dis_mg$model$pty_dis))))

pr_unity_int <- summary(prediction(pr_uty_dis_id_int_mg,
                               vcov = vcov_cluster(pr_uty_dis_id_int, cluster = ~agency),
                               at = list(ideo_rating = unique(pr_uty_dis_id_int_mg$model$ideo_rating),
                                         pty_dis = c(-1.25, 0, 1.5))))

pr_vv_int <- summary(prediction(pr_vv_dis_id_int_mg,
                               vcov = vcov_cluster(pr_vv_dis_id_int, cluster = ~agency),
                               at = list(ideo_rating = unique(pr_vv_dis_id_int_mg$model$ideo_rating),
                                         pty_dis = c(-1.25, 0, 1.5))))

#### Figure A.8 ####

pl_unity <- ggplot(pr_unity, aes(`at(pty_dis)`, Prediction)) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(x = `at(pty_dis)`,
                  ymin = lower, ymax = upper),
                  alpha = 0.25) + 
  labs(title = "Party Unity Votes",
       x = "Partisan Disagreement",
       y = "Pred. Prob. of Party Unity Vote") +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1),
                     limits = c(-0.05, 1.1)) +
  theme_bw() +
  theme(legend.position = "none")

pl_unity_int <-  ggplot(pr_unity_int, aes(`at(ideo_rating)`, Prediction,
                                          shape = as.factor(`at(pty_dis)`),
                                          color = as.factor(`at(pty_dis)`),
                                          fill = as.factor(`at(pty_dis)`))) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(x = `at(ideo_rating)`,
                  ymin = lower, ymax = upper),
              alpha = 0.25, show.legend = FALSE) + 
  labs(x = "Agency Ideology",
       y = "Pred. Prob. of Party Unity Vote",
       shape = "Partisan Disagreement",
       color = "Partisan Disagreement",
       fill = "Partisan Disagreement") +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1),
                     limits = c(-0.2, 1.2)) +
  scale_shape_manual(values = 15:18,
                     labels = c("Low", "Moderate", "High")) +
  scale_color_manual(values = c("gray57", "orange", "purple"),
                     labels = c("Low", "Moderate", "High")) +
  scale_fill_manual(values = c("gray57", "orange", "purple"),
                    labels = c("Low", "Moderate", "High")) +
  theme_bw() +
  theme(legend.position = "none")

pl_vv <- ggplot(pr_vv, aes(`at(pty_dis)`, Prediction)) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(x = `at(pty_dis)`,
                  ymin = lower, ymax = upper),
              alpha = 0.25) + 
  labs(title = "Voice Votes",
       x = "Partisan Disagreement",
       y = "Pred. Prob. of Voice Vote") +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1),
                     limits = c(-0.05, 1.1)) +
  theme_bw() +
  theme(legend.position = "none")

pl_vv_int <- ggplot(pr_vv_int, aes(`at(ideo_rating)`, Prediction,
                                          shape = as.factor(`at(pty_dis)`),
                                          color = as.factor(`at(pty_dis)`),
                                          fill = as.factor(`at(pty_dis)`))) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(x = `at(ideo_rating)`,
                  ymin = lower, ymax = upper),
              alpha = 0.25, show.legend = FALSE) + 
  labs(x = "Agency Ideology",
       y = "Pred. Prob of Voice Vote",
       shape = "Partisan Disagreement",
       color = "Partisan Disagreement",
       fill =  "Partisan Disagreement") +
  scale_y_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1),
                     limits = c(-0.2, 1.2)) +
  scale_shape_manual(values = 15:18,
                     labels = c("Low", "Moderate", "High")) +
  scale_color_manual(values = c("gray57", "orange", "purple"),
                     labels = c("Low", "Moderate", "High")) +
  scale_fill_manual(values = c("gray57", "orange", "purple"),
                    labels = c("Low", "Moderate", "High")) + 
  theme_bw()

(pl_unity + pl_vv) / (pl_unity_int + pl_vv_int) +
  plot_annotation(caption = "Shaded area gives the 95% confidence interval.") +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

# Clean up
rm(list = ls()[!str_detect(ls(), "dis_data_wt")])

#### Section A.8.1 ####

setFixest_dict(pty_dis = "Partisan Disagreement",
               ideo_rating = "Agency Ideology",
               `abs(ideo_rating)` = "|Agency Ideology|",
               `I(ideo_rating^2)` = "Agency Ideology Squared")

# Fit the models
lm_rcl_1 <- feols(pty_dis ~ ideo_rating,
                  data = dis_data_wt, vcov = "hetero")

lm_rcl_2 <- feols(pty_dis ~ abs(ideo_rating),
                    data = dis_data_wt, vcov = "hetero")

lm_rcl_3 <- feols(pty_dis ~ ideo_rating + I(ideo_rating^2),
                    data = dis_data_wt, vcov = "hetero")

# Standard errors of the regressions

round(sqrt(lm_rcl_1$sigma2), digits = 2)
round(sqrt(lm_rcl_3$sigma2), digits = 2)
round(sqrt(lm_rcl_2$sigma2), digits = 2)

#### Table A.6 ####

etable(lm_rcl_1, lm_rcl_3, lm_rcl_2,
       tex = TRUE, digits = "r2",
       digits.stats = "r2",
       order = "!Constant",
       style.tex = style.tex(stats.title = "\\midrule"),
       title = "Predicting Partisan Disagreement with Agency Ideology")



#### Figure A.9 ####

# Create lm objects for prediction()
lm_2 <- lm(pty_dis ~ abs_ideo,
           data = dis_data_wt %>% mutate(abs_ideo = abs(ideo_rating)))

lm_3 <- lm(pty_dis ~ ideo_rating + I(ideo_rating^2),
           data = dis_data_wt)


# Note: feols() uses HC1
pr_abs <- summary(prediction(lm_2,
                             vcov = sandwich::vcovHC(lm_2, type = "HC1"),
                             at = list(abs_ideo = unique(lm_2$model$abs_ideo))))

pr_quad <- summary(prediction(lm_3,
                              vcov = sandwich::vcovHC(lm_3, type = "HC1"),
                              at = list(ideo_rating = unique(lm_3$model$ideo_rating))))


# Create the figure

quad <- ggplot() +
  geom_point(data = dis_data_wt,
             mapping = aes(ideo_rating, pty_dis,
                           shape = ideo_cat, color = ideo_cat),
             size = 2, alpha = 1) +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  geom_line(data = pr_quad, mapping = aes(`at(ideo_rating)`, Prediction)) +
  geom_ribbon(data = pr_quad, mapping = aes(x = `at(ideo_rating)`,
                                       ymin = lower, ymax = upper),
              alpha = 0.25) +
  labs(x = "Agency Ideology",
       y = "Partisan Disagreement",
       color = "Agency Ideology",
       shape = "Agency Ideology") +
  theme_bw()

abs <- ggplot() +
  geom_point(data = dis_data_wt,
             mapping = aes(abs(ideo_rating), pty_dis,
                           shape = ideo_cat, color = ideo_cat),
             size = 2, alpha = 1) +
  scale_shape_manual(values = 15:18) +
  scale_color_manual(values = c("firebrick2", "gray57", "dodgerblue", "black")) +
  geom_line(data = pr_abs, mapping = aes(`at(abs_ideo)`, Prediction)) +
  geom_ribbon(data = pr_abs, mapping = aes(x = `at(abs_ideo)`,
                                       ymin = lower, ymax = upper),
              alpha = 0.25) +
  labs(x = "|Agency Ideology|",
       y = "Partisan Disagreement",
       color = "Agency Ideology",
       shape = "Agency Ideology") +
  theme_bw()

(quad + abs) +
  plot_annotation(caption = "The solid lines are predicted partisan disagreement from the quadratic and absolute-value models in Table A.6.\n The shaded region gives the 95% confidence interval.") +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")


